High Performance Cholesky Factorization via Blocking and Recursion That Uses Minimal Storage

We present a high performance Cholesky factorization algorithm , called BPC for Blocked Packed Cholesky, which performs better or equivalent to the LAPACK DPOTRF subroutine, but with about the same memory requirements as the LAPACK DPPTRF subroutine, which runs at level 2 BLAS speed. Algorithm BPC only calls DGEMM and level 3 kernel routines. It combines a recursive algorithm with blocking and a recursive packed data format. A full analysis of overcoming the non-linear addressing overhead imposed by recursion is given and discussed. Finally, since BPC uses GEMM to a great extent, we easily get a considerable amount of SMP parallelism from an SMP GEMM. In 1], a new recursive packed format for Cholesky factorization was described. It is a variant of the triangular format described in 5], and suggested in 4]. For a reference to recursion in linear algebra, see 4, 5]. The advantage of using packed formats instead of full format is the memory savings possible when working with large matrices. The disadvantage is that the use of high-performance standard library routines such as the matrix multiply and add operation, GEMM, is inhibited. We combine blocking with recursion to produce a practical implementation of a new algorithm called BPC, for blocked packed Cholesky 3, pages 142-147]. BPC rst transforms standard packed lower format to packed recursive lower row format. The during execution algorithm BPC only calls standard DGEMM and level 3 kernel routines. This method has the beneet of being transparent to the user. No extra assumptions about the input matrix needs be made. This is important since an algorithm like Cholesky using a new data format would not work on previous codes. Algorithm BPC outperforms level 3 LAPACK 2] routine DPORTF on three IBM platforms, the POWER3, the POWER2, and the POWERPC 604e, even when